Eggshell_review_project
Lucas Legendre
25/04/2022
Compiled under R version 4.2.0 (2022-04-22)
WARNING: edit the working directory to your preferred folder.
This document details all analyses performed in R for the
study:
Legendre, L. J., S. Choi, and J. A. Clarke. The use of inconsistent
terminology for reptile eggshell traits affects the outcome of
evolutionary analyses. Journal of Anatomy.
For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.
Loading packages
library(ape)
library(castor)
library(evobiR)
library(ggtree)
library(phytools)
library(RColorBrewer)Data and tree
- Tree
tree<-read.nexus("treewhole_newversion.trees.nex")
plotTree(tree, fsize=0.5,lwd=1,type="fan",ftype="i")
- Data
data<-read.table("Datawhole_newproject.txt", header=T)
setdiff(tree$tip.label, data$Taxon) # taxa and data match## character(0)rownames(data)<-data$Taxon
datanew<-ReorderData(tree, data)Ancestral state reconstruction: discrete trait (soft/semi-rigid/hard)
- Extract data vectors for each coding
thisstudy<-datanew$Type_2021; names(thisstudy)<-rownames(datanew)
norell<-datanew$Type_Norell; names(norell)<-rownames(datanew)
legendre<-datanew$Type_Legendre; names(legendre)<-rownames(datanew)
# Color palette
cols<-setNames(c("royalblue","green3","red3"),c("Soft","Semi-rigid","Hard"))
# Visualize tree with node numbers
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
Using character coding defined in this study (new scoring)
- Perform SIMMAP using
phytools– 1000 iterations, using AIC to select best model (code modified from Liam Revell, see here)
x<-thisstudy
aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
d.aic<-aic-min(aic)
exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.001035456 0.000517728 0.000517728
## Semi-rigid 0.000517728 -0.001035456 0.000517728
## Soft 0.000517728 0.000517728 -0.001035456
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010899924 0.0004565271 0.0006334653
## Semi-rigid 0.0004565271 -0.0007161171 0.0002595901
## Soft 0.0006334653 0.0002595901 -0.0008930554
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.0000000000 0.0000000000
## Semi-rigid 0.0072145631 -0.0109014209 0.0036868578
## Soft 0.0004663626 0.0002355611 -0.0007019237
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.logL## ER SYM ARD
## -56.68015 -56.39465 -45.61188AIC<-mapply(aic,logL,c(1,3,6))
AIC## ER SYM ARD
## 115.3603 118.7893 103.2238AIC.W<-aic.w(AIC)
AIC.W## ER SYM ARD
## 0.0023088700 0.0004157184 0.9972754116nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 2 1 997o1<-make.simmap(tree,x,model="ARD",nsim=nsim[1])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.0000000000 0.0000000000
## Semi-rigid 0.0072145631 -0.0109014209 0.0036868578
## Soft 0.0004663626 0.0002355611 -0.0007019237
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.o2<-make.simmap(tree,x,model="ER",nsim=nsim[2])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.001035456 0.000517728 0.000517728
## Semi-rigid 0.000517728 -0.001035456 0.000517728
## Soft 0.000517728 0.000517728 -0.001035456
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treesstudy<-c(o1,o2)
objstudy<-describe.simmap(treesstudy)
plot(objstudy,type="fan",fsize=0.01,lwd=1,ftype="i", colors=cols,ylim=c(-2,Ntip(tree)),offset=20, part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Here, a semi-rigid eggshell is the ancestral state for almost all major reptile clades – reptiles, lepidosaurs, squamates, turtles, archelosaurs, archosaurs, dinosaurs,ornithischians, and saurischians. However, there are only 6 taxa with semi-rigid eggshells out of 208!
=> Some of these taxa may have a excessive influence on this result – probably the two sauropodomorphs with that eggshell type, since they are the closer in age to the root and to the aforementioned subclades.
Using character coding from Norell et al. (2020) – ratio scoring
x<-norell
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013917968 0.0006958984 0.0006958984
## Semi-rigid 0.0006958984 -0.0013917968 0.0006958984
## Soft 0.0006958984 0.0006958984 -0.0013917968
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013792703 0.0004831551 0.0008961152
## Semi-rigid 0.0004831551 -0.0010895186 0.0006063635
## Soft 0.0008961152 0.0006063635 -0.0015024787
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.000000000 0.000000000 0.000000000
## Semi-rigid 0.008810469 -0.014180320 0.005369851
## Soft 0.001156338 0.001443489 -0.002599828
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.logL## ER SYM ARD
## -70.17246 -69.66276 -57.83996AIC<-mapply(aic,logL,c(1,3,6))
AIC## ER SYM ARD
## 142.3449 145.3255 127.6799AIC.W<-aic.w(AIC)
AIC.W## ER SYM ARD
## 0.0006534088 0.0001472160 0.9991993752nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 1 0 999n1<-make.simmap(tree,x,model="ARD",nsim=nsim[1])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.000000000 0.000000000 0.000000000
## Semi-rigid 0.008810469 -0.014180320 0.005369851
## Soft 0.001156338 0.001443489 -0.002599828
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.n2<-make.simmap(tree,x,model="ER",nsim=nsim[2])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013917968 0.0006958984 0.0006958984
## Semi-rigid 0.0006958984 -0.0013917968 0.0006958984
## Soft 0.0006958984 0.0006958984 -0.0013917968
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treesnorell<-c(o1,o2)
objnorell<-describe.simmap(treesnorell)
plot(objnorell,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Completely different result – let us check which taxa are different from the previous ASR…
datanew[c(which(datanew$Type_2021!=datanew$Type_Norell)),]# Only 7 taxa are different – two non-avian dinosaurs and five turtles
# Changing character state for the two non-avian dinosaurs
x[c(21,22)]<-"Semi-rigid"
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0014141228 0.0007070614 0.0007070614
## Semi-rigid 0.0007070614 -0.0014141228 0.0007070614
## Soft 0.0007070614 0.0007070614 -0.0014141228
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013866412 0.0005855999 0.0008010413
## Semi-rigid 0.0005855999 -0.0013144727 0.0007288728
## Soft 0.0008010413 0.0007288728 -0.0015299141
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.0000000000 0.000000000
## Semi-rigid 0.0081716323 -0.0131662651 0.004994633
## Soft 0.0005407866 0.0007809644 -0.001321751
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 0 0 1000treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.0000000000 0.000000000
## Semi-rigid 0.0081716323 -0.0131662651 0.004994633
## Soft 0.0005407866 0.0007809644 -0.001321751
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.objnorell<-describe.simmap(treesnorell)
plot(objnorell,fsize=0.5,lwd=1,ftype="i",colors=cols,ylim=c(-2,Ntip(tree)),part=0.95)
add.simmap.legend(colors=cols,prompt=FALSE,x=20,y=180)
All major clades are now recovered as semi-rigid again => strong influence of the two sauropodomorphs.
Using character coding from Legendre et al. (2020) – shell unit scoring
x<-legendre
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007487292 0.0003743646 0.0003743646
## Semi-rigid 0.0003743646 -0.0007487292 0.0003743646
## Soft 0.0003743646 0.0003743646 -0.0007487292
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006125127 0.0000000000 0.0006125127
## Semi-rigid 0.0000000000 -0.0005646183 0.0005646183
## Soft 0.0006125127 0.0005646183 -0.0011771310
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004334394 0.000000000 0.0004334394
## Semi-rigid 0.0034884449 -0.007010172 0.0035217272
## Soft 0.0003883791 0.000000000 -0.0003883791
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.logL## ER SYM ARD
## -45.68062 -43.28735 -41.37627AIC<-mapply(aic,logL,c(1,3,6))
AIC## ER SYM ARD
## 93.36124 92.57469 94.75254AIC.W<-aic.w(AIC)
AIC.W## ER SYM ARD
## 0.3355056 0.4971605 0.1673339nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 336 497 167l1<-make.simmap(tree,x,model="ER",nsim=nsim[1])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007487292 0.0003743646 0.0003743646
## Semi-rigid 0.0003743646 -0.0007487292 0.0003743646
## Soft 0.0003743646 0.0003743646 -0.0007487292
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.l2<-make.simmap(tree,x,model="SYM",nsim=nsim[2])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006125127 0.0000000000 0.0006125127
## Semi-rigid 0.0000000000 -0.0005646183 0.0005646183
## Soft 0.0006125127 0.0005646183 -0.0011771310
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.l3<-make.simmap(tree,x,model="ARD",nsim=nsim[3])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004334394 0.000000000 0.0004334394
## Semi-rigid 0.0034884449 -0.007010172 0.0035217272
## Soft 0.0003883791 0.000000000 -0.0003883791
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treeslegendre<-c(l1,l2,l3)
objlegendre<-describe.simmap(treeslegendre)
plot(objlegendre,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Reptiles, lepidosaurs, and squamates are recovered as ancestrally soft-shelled, while turtles, archelosaurs, archosaurs, and all dinosaur clades are recovered as hard-shelled.
The pattern seems to follow whichever character state is coded for the two sauropodomorphs Lufengosaurus and Massospondylus.
To test this hypothesis, we must remove all other taxa (n = 7) that are not coded the same in all three studies, and see if the pattern is still present.
- Remove all taxa with differences in coding, except the two sauropodomorphs
remove<-rownames(datanew[c(which(datanew$Type_Legendre!=datanew$Type_Norell)),][c(3:9),])
x<-thisstudy[!names(thisstudy)%in%remove]
treenew<-drop.tip(tree, setdiff(tree$tip.label, names(x)))
# Visualize tree with node numbers
ggtree(treenew) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
Replicate analysis on the new tree 3 times, with identical coding for all taxa except the two sauropodomorphs
- Coded as semi-rigid
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008797485 0.0004398742 0.0004398742
## Semi-rigid 0.0004398742 -0.0008797485 0.0004398742
## Soft 0.0004398742 0.0004398742 -0.0008797485
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -7.696326e-04 6.160148e-05 0.0007080311
## Semi-rigid 6.160148e-05 -7.755548e-04 0.0007139533
## Soft 7.080311e-04 7.139533e-04 -0.0014219844
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0002017335 0.0000000000 0.0002017335
## Semi-rigid 0.0037018375 -0.0037018375 0.0000000000
## Soft 0.0020196421 0.0008852088 -0.0029048508
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 388 195 417s1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008797485 0.0004398742 0.0004398742
## Semi-rigid 0.0004398742 -0.0008797485 0.0004398742
## Soft 0.0004398742 0.0004398742 -0.0008797485
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.s2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -7.696326e-04 6.160148e-05 0.0007080311
## Semi-rigid 6.160148e-05 -7.755548e-04 0.0007139533
## Soft 7.080311e-04 7.139533e-04 -0.0014219844
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.s3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0002017335 0.0000000000 0.0002017335
## Semi-rigid 0.0037018375 -0.0037018375 0.0000000000
## Soft 0.0020196421 0.0008852088 -0.0029048508
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treenewsr<-c(s1,s2,s3)
objsr<-describe.simmap(treenewsr)
plot(objsr,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Most inner clades are now recovered as soft-shelled – most likely an artefact due to the position of Mussaurus – closer to the other inner nodes than any other taxon, since it is the oldest specimen in the sample.
- Coded as soft
x[c(21,22)]<-"Soft"
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008695627 0.0004347813 0.0004347813
## Semi-rigid 0.0004347813 -0.0008695627 0.0004347813
## Soft 0.0004347813 0.0004347813 -0.0008695627
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007415283 0.0000000000 0.0007415283
## Semi-rigid 0.0000000000 -0.0005402947 0.0005402947
## Soft 0.0007415283 0.0005402947 -0.0012818231
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001997071 0.0000000000 0.0001997071
## Semi-rigid 0.0023041591 -0.0023041591 0.0000000000
## Soft 0.0021232656 0.0005497209 -0.0026729866
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 113 284 603so1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008695627 0.0004347813 0.0004347813
## Semi-rigid 0.0004347813 -0.0008695627 0.0004347813
## Soft 0.0004347813 0.0004347813 -0.0008695627
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.so2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007415283 0.0000000000 0.0007415283
## Semi-rigid 0.0000000000 -0.0005402947 0.0005402947
## Soft 0.0007415283 0.0005402947 -0.0012818231
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.so3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001997071 0.0000000000 0.0001997071
## Semi-rigid 0.0023041591 -0.0023041591 0.0000000000
## Soft 0.0021232656 0.0005497209 -0.0026729866
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treenews<-c(so1, so2, so3)
objs<-describe.simmap(treenews)
plot(objs,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Same result as with semi-rigid coding, unsurprisingly.
- Coded as hard
x[c(21,22)]<-"Hard"
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007710002 0.0003855001 0.0003855001
## Semi-rigid 0.0003855001 -0.0007710002 0.0003855001
## Soft 0.0003855001 0.0003855001 -0.0007710002
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006317588 0.000000000 0.0006317588
## Semi-rigid 0.0000000000 -0.000568647 0.0005686470
## Soft 0.0006317588 0.000568647 -0.0012004058
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004499743 0.000000000 0.0004499743
## Semi-rigid 0.0034859171 -0.007000114 0.0035141964
## Soft 0.0003928020 0.000000000 -0.0003928020
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 347 497 156h1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007710002 0.0003855001 0.0003855001
## Semi-rigid 0.0003855001 -0.0007710002 0.0003855001
## Soft 0.0003855001 0.0003855001 -0.0007710002
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.h2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006317588 0.000000000 0.0006317588
## Semi-rigid 0.0000000000 -0.000568647 0.0005686470
## Soft 0.0006317588 0.000568647 -0.0012004058
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.h3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0004499743 0.000000000 0.0004499743
## Semi-rigid 0.0034859171 -0.007000114 0.0035141964
## Soft 0.0003928020 0.000000000 -0.0003928020
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treenewh<-c(h1, h2, h3)
objh<-describe.simmap(treenewh)
plot(objh,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Ancestral state for dinosaurs and archosaurs, and all major less inclusive clades except pterosaurs, is hard-shelled, as in the Legendre et al. (2020) coding.
To check how strongly are these results influenced by branch length
information, we can replicate these analyses using maximum parsimony,
which does not consider branch length information, using
castor.
- Coding from this study
studyN<-thisstudy
studyN[studyN=="Soft"]<-1; studyN[studyN=="Semi-rigid"]<-2; studyN[studyN=="Hard"]<-3
studyN<-as.numeric(studyN); names(studyN)<-names(thisstudy)
MPS<-asr_max_parsimony(tree, studyN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPS$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(studyN,sort(unique(studyN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)
# ancestral states
MPS2<-MPS$ancestral_likelihoods
rownames(MPS2)<-c(209:374); colnames(MPS2)<-c("soft","semi-rigid","hard")- Coding from Norell et al. (2020)
norellN<-norell
norellN[norellN=="Soft"]<-1; norellN[norellN=="Semi-rigid"]<-2; norellN[norellN=="Hard"]<-3
norellN<-as.numeric(norellN); names(norellN)<-names(norell)
MPN<-asr_max_parsimony(tree, norellN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPN$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(norellN,sort(unique(norellN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=180,fsize=0.8)
# ancestral states
MPN2<-MPN$ancestral_likelihoods
rownames(MPN2)<-c(209:374); colnames(MPN2)<-c("soft","semi-rigid","hard")- Coding from Legendre et al. (2020)
legendreN<-legendre
legendreN[legendreN=="Soft"]<-1; legendreN[legendreN=="Semi-rigid"]<-2; legendreN[legendreN=="Hard"]<-3
legendreN<-as.numeric(legendreN); names(legendreN)<-names(legendre)
MPL<-asr_max_parsimony(tree, legendreN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPL$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(legendreN,sort(unique(legendreN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)
# ancestral states
MPL2<-MPL$ancestral_likelihoods
rownames(MPL2)<-c(209:374); colnames(MPL2)<-c("soft","semi-rigid","hard")- With the reduced tree, changing only the two sauropodomorphs
xN<-x
xN[xN=="Soft"]<-1; xN[xN=="Semi-rigid"]<-2; xN[xN=="Hard"]<-3
xN<-as.numeric(xN); names(xN)<-names(x)
# Coded as soft
xN[c(21,22)]<-1
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")
# Coded as semi-rigid
xN[c(21,22)]<-2
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")
# Coded as hard
xN[c(21,22)]<-3
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")Ancestral state reconstruction – continuous character (eggshell thickness)
- Data and tree
treedi<-multi2di(tree, random=FALSE)
# To visualize node numbers on the tree
ggtree(treedi) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
AET<-log1p(datanew$Eggshell_thickness); names(AET)<-rownames(datanew)
RET<-log1p(datanew$Eggshell_thickness/datanew$Egg_mass); names(RET)<-rownames(datanew)- Phylogenetic signal (
contMapassumes a Brownian motion model)
phylosig(treedi, AET, method="lambda", test=TRUE)##
## Phylogenetic signal lambda : 0.999934
## logL(lambda) : -316.25
## LR(lambda=0) : 255.021
## P-value (based on LR test) : 2.08867e-57phylosig(treedi, RET, method="lambda", test=TRUE)##
## Phylogenetic signal lambda : 0.999934
## logL(lambda) : -231.389
## LR(lambda=0) : 206.653
## P-value (based on LR test) : 7.38001e-47Very strong signal for both absolute and relative eggshell thickness.
- ASR for absolute eggshell thickness
phyloplotAET<-contMap(treedi, AET, plot=F)
plot(setMap(phyloplotAET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
prompt=FALSE,x=-300,y=280)
# Ancestral states
expm1(fastAnc(treedi, AET))## Ancestral character estimates using fastAnc:
## 209 210 211 212 213 214
## 30.237176 30.237176 27.637279 20.146957 37.147079 57.434822
## 215 216 217 218 219 220
## 59.490746 38.962770 37.947764 62.376276 63.738529 60.890954
## 221 222 223 224 225 226
## 18.797940 17.535632 17.058187 21.413045 14.494621 17.213248
## 227 228 229 230 231 232
## 16.971098 16.477399 14.535149 16.143759 8.220722 8.894198
## 233 234 235 236 237 238
## 17.611084 17.663373 11.373368 11.394027 11.094212 7.685149
## 239 240 241 242 243 244
## 18.197125 27.804160 11.643952 6.395467 5.262258 11.872530
## 245 246 247 248 249 250
## 10.876488 12.263387 12.123808 14.013434 13.895900 14.226852
## 251 252 253 254 255 256
## 16.558085 31.172635 40.914364 133.738678 544.568527 573.610635
## 257 258 259 260 261 262
## 610.795533 607.951060 498.450624 137.074853 142.357447 143.074080
## 263 264 265 266 267 268
## 139.614328 143.141126 61.052662 155.004338 178.337602 198.598257
## 269 270 271 272 273 274
## 196.829154 167.286930 145.639082 217.263595 259.777176 324.061138
## 275 276 277 278 279 280
## 38.519274 253.555625 279.270796 414.658851 624.034336 367.028433
## 281 282 283 284 285 286
## 460.030533 450.070207 468.619867 478.051897 33.579247 2.962142
## 287 288 289 290 291 292
## 0.014681 4.087266 11.285567 3.380069 41.916492 152.227600
## 293 294 295 296 297 298
## 152.227600 152.227600 1426.916211 1426.916211 152.227600 1818.852274
## 299 300 301 302 303 304
## 1818.852274 1455.301611 1501.790647 31.092007 14.703935 72.646286
## 305 306 307 308 309 310
## 1.813370 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348
## 311 312 313 314 315 316
## 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348
## 317 318 319 320 321 322
## 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348
## 323 324 325 326 327 328
## 1725.926348 321.949973 479.318439 856.292119 856.292119 856.292119
## 329 330 331 332 333 334
## 533.581298 535.746398 535.746398 557.800038 1286.491850 567.795787
## 335 336 337 338 339 340
## 606.755787 606.755787 606.755787 606.755787 1391.768701 1101.675921
## 341 342 343 344 345 346
## 540.855621 540.855621 540.855621 656.650226 418.648969 418.648969
## 347 348 349 350 351 352
## 418.648969 400.036092 400.036092 400.036092 400.036092 383.759407
## 353 354 355 356 357 358
## 393.341060 482.609908 482.609908 482.609908 482.609908 482.609908
## 359 360 361 362 363 364
## 487.383928 482.609908 672.440036 688.378914 743.268097 749.865095
## 365 366 367 368 369 370
## 642.582251 1304.089170 319.366201 778.769338 908.526782 820.791309
## 371 372 373 374 375 376
## 423.102753 255.042210 300.656066 343.613100 298.832008 290.376178
## 377 378 379 380 381 382
## 222.124964 217.737359 177.855640 166.738135 161.028661 229.680422
## 383 384 385 386 387 388
## 234.234604 236.215710 233.547118 232.188904 231.149423 196.464412
## 389 390 391 392 393 394
## 216.848165 210.911249 208.760967 210.718962 206.849992 201.750368
## 395 396 397 398 399 400
## 251.215833 289.694249 316.988394 253.496044 257.536715 243.363765
## 401 402 403 404 405 406
## 151.833016 155.871496 164.021661 144.782394 128.404579 106.959834
## 407 408 409 410 411 412
## 100.296809 118.839823 109.667780 73.949894 103.667666 99.285850
## 413 414 415
## 98.490499 103.797577 101.497881Ancestral eggshell thickness seems to be intermediate for all major nodes, including archosaurs and dinosaurs. However, many extant lepidosaurs are also recovered as intermediate – likely due to the extremely low values in pterosaurs, reported as lacking a calcareous layer entirely, which shift any thicker eggshell towards the middle of the spectrum.
Furthermore, the pattern strongly follows that of body mass, as expected (e.g. Stein et al., 2019; Legendre and Clarke, 2021).
- ASR for relative eggshell thickness
phyloplotRET<-contMap(treedi, RET, plot=F)
plot(setMap(phyloplotRET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
prompt=FALSE,x=-300,y=280)
# Ancestral states
expm1(fastAnc(treedi, RET))## Ancestral character estimates using fastAnc:
## 209 210 211 212 213 214 215 216
## 3.463807 3.463807 4.072394 7.776137 20.371286 37.945868 39.830943 25.691073
## 217 218 219 220 221 222 223 224
## 39.657850 42.685185 44.863801 53.013755 7.517838 7.370333 7.179543 6.555069
## 225 226 227 228 229 230 231 232
## 9.463869 7.389222 7.661106 9.194875 11.651921 13.025134 12.735303 13.646880
## 233 234 235 236 237 238 239 240
## 16.532983 20.147667 21.875358 24.370764 23.821319 14.769232 1.647382 1.003291
## 241 242 243 244 245 246 247 248
## 1.610748 0.936602 0.860325 1.828704 1.596699 2.012611 3.474841 1.683005
## 249 250 251 252 253 254 255 256
## 1.180153 1.164833 2.085560 3.133937 2.664801 6.649956 10.831692 10.811365
## 257 258 259 260 261 262 263 264
## 10.536241 10.641737 9.256815 6.921701 8.109834 10.716757 10.797700 12.176394
## 265 266 267 268 269 270 271 272
## 2.043793 11.104335 10.749829 10.936284 10.273487 9.993487 8.767791 10.221733
## 273 274 275 276 277 278 279 280
## 7.286859 7.960865 2.192523 4.759820 4.819961 5.480688 6.780228 4.590287
## 281 282 283 284 285 286 287 288
## 5.810148 6.065964 5.670210 5.486058 1.755410 0.464691 0.004048 0.554308
## 289 290 291 292 293 294 295 296
## 1.153786 0.409454 1.322077 1.359789 1.359789 1.359789 5.337305 5.337305
## 297 298 299 300 301 302 303 304
## 1.359789 0.659302 0.659302 1.233747 1.279996 0.923542 0.437393 0.389384
## 305 306 307 308 309 310 311 312
## 0.131437 0.730699 0.730699 0.730699 0.730699 0.730699 0.730699 0.730699
## 313 314 315 316 317 318 319 320
## 0.730699 0.730699 0.730699 0.730699 0.730699 0.730699 0.730699 0.730699
## 321 322 323 324 325 326 327 328
## 0.730699 0.730699 0.730699 2.798346 4.269776 4.207026 4.207026 4.207026
## 329 330 331 332 333 334 335 336
## 5.096371 5.317186 5.317186 5.249324 2.804855 5.265575 4.919952 4.919952
## 337 338 339 340 341 342 343 344
## 4.919952 4.919952 2.041476 3.280823 5.648622 5.648622 5.648622 4.639335
## 345 346 347 348 349 350 351 352
## 8.397269 8.397269 8.397269 8.721975 8.721975 8.721975 8.721975 8.944506
## 353 354 355 356 357 358 359 360
## 8.806187 8.062933 8.062933 8.062933 8.062933 8.062933 8.107327 8.062933
## 361 362 363 364 365 366 367 368
## 3.124942 3.017798 2.342437 2.259183 2.383122 0.661661 4.634284 2.112206
## 369 370 371 372 373 374 375 376
## 1.859723 1.878858 1.176578 10.834138 9.311660 7.858083 9.853769 10.770714
## 377 378 379 380 381 382 383 384
## 11.860656 11.582529 13.304265 13.852859 14.336956 10.818602 10.468547 10.019546
## 385 386 387 388 389 390 391 392
## 9.677966 9.355086 9.031513 13.464317 11.100694 10.828751 10.901636 10.572284
## 393 394 395 396 397 398 399 400
## 10.469217 10.385413 9.574631 8.311346 7.455384 9.415512 8.811797 9.272837
## 401 402 403 404 405 406 407 408
## 19.516988 20.275505 19.452318 20.563944 23.256399 28.370771 30.397447 25.081937
## 409 410 411 412 413 414 415
## 27.482967 41.577914 29.373696 30.524621 32.153419 31.108670 34.366359Much lower values for archosaurs and dinosaurs – eggshell thickness seems to have become thinner for a given egg mass at the base of Ornithodira, and later increased in theropods. However, this pattern is dependent on a very small sample of pterosaurs, ornithischians, and sauropods; it is likely that the addition of new specimens attributed to any of these clades would considerably change that pattern.
We can see a strong increase in geckos and in eufalconimorphs – the latter having already been identified in Legendre and Clarke (2021).
References
- Legendre, L.J., Clarke, J.A., 2021. Shifts in eggshell thickness are related to changes in locomotor ecology in dinosaurs. Evolution 75, 1415–1430.
- Legendre, L.J., Rubilar-Rogers, D., Musser, G.M., Davis, S.N., Otero, R.A., Vargas, A.O., Clarke, J.A., 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature 583, 411–414.
- Norell, M.A., Wiemann, J., Fabbri, M., Yu, C., Marsicano, C.A., Moore-Nall, A., Varricchio, D.J., Pol, D., Zelenitsky, D.K., 2020. The first dinosaur egg was soft. Nature 583, 406–410.
- Stein, K., Prondvai, E., Huang, T., Baele, J.-M., Sander, P.M., Reisz, R., 2019. Structure and evolutionary implications of the earliest (Sinemurian, Early Jurassic) dinosaur eggs and eggshells. Scientific Reports 9, 4424.